suppressPackageStartupMessages({
require(tidyr)
require(lme4) # for running models
require(sjPlot) # for visualizing model interactions
require(sjmisc)
require(sjlabelled)
require(tidyverse)
require(viridis)
require(lmerTest)
})

Climate data

Population_climateNA <- read.csv("C:/Dropbox/UVM/Research/Code/ClimateData/ClimateNA/Family/PCA_sprucePops.txt",
                                 header=T, sep="\t")
colnames(Population_climateNA)[colnames(Population_climateNA)=="PC1_sel"] <- "PC1"

Family_climateNA <- read.csv("C:/Dropbox/UVM/Research/Code/ClimateData/ClimateNA/Family/PCA_spruceFams.txt",
                                 header=T, sep="\t")
colnames(Family_climateNA)[colnames(Family_climateNA)=="PC1_sel"] <- "PC1"

Trait data

budset_2019 <- read.csv("./trait_data/BudSet_2019.csv", stringsAsFactors = T)
budset_2020 <- read.csv("./trait_data/BudSet_2020.csv", stringsAsFactors = T)
budbreak_2020 <- read.csv("./trait_data/BudBreak_2020_cGDD.csv", stringsAsFactors = T)
growth_2019 <- read.csv("./trait_data/Growth_2019.csv", stringsAsFactors = T)
growth_2020 <- read.csv("./trait_data/Growth_2020.csv", stringsAsFactors = T)
Plasticity <- read.csv("./trait_data/Plasticity_TD.csv", stringsAsFactors = T)
  • Reorder levels
# Garden  
budset_2019$Garden <-  factor(budset_2019$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
budset_2020$Garden <-  factor(budset_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
budbreak_2020$Garden <-  factor(budbreak_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
growth_2019$Garden <-  factor(growth_2019$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
growth_2020$Garden <-  factor(growth_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))

# Region
budset_2019$Region <-  factor(budset_2019$Region,levels = c("Core", "Margin", "Edge"))
budset_2020$Region <-  factor(budset_2020$Region,levels = c("Core", "Margin", "Edge"))
budbreak_2020$Region <-  factor(budbreak_2020$Region,levels = c("Core", "Margin", "Edge"))
growth_2019$Region <-  factor(growth_2019$Region,levels = c("Core", "Margin", "Edge"))
growth_2020$Region <-  factor(growth_2020$Region,levels = c("Core", "Margin", "Edge"))

Models

Growth 2019

## working model for growth - ****
Growth19_mod1 <- lmer(data = growth_2019, na.action = na.omit,
               Growth ~ Garden * Region +
               (1|Population) + (1|mBed) +
               (Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
Growth19_mod1a <- lmer(data = growth_2019, na.action = na.omit,
               Growth ~ Garden * Region +
               (1|Population) + (1|mBed) + (1|Family))

Growth19_mod1b <- lmer(data = growth_2019, na.action = na.omit,
               Growth ~ Garden * Region +
               (1|Population) + (1|mBed))

# model 1 anova 
anova(Growth19_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Garden          75.68   37.84     2  12.32  1.0818  0.369101    
## Region        2991.90 1495.95     2  64.21 42.7664 1.562e-12 ***
## Garden:Region  532.33  133.08     4 722.10  3.8046  0.004536 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(Growth19_mod1, Growth19_mod1a) 
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## Growth19_mod1a:     Family)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)   
## Growth19_mod1a   13 30358 30442 -15166    30332                        
## Growth19_mod1    18 30351 30467 -15158    30315 16.499  5   0.005555 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(Growth19_mod1a, Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## Growth19_mod1a:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## Growth19_mod1b   12 30469 30547 -15223    30445                         
## Growth19_mod1a   13 30358 30442 -15166    30332 113.39  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(Growth19_mod1, Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## Growth19_mod1b   12 30469 30547 -15223    30445                         
## Growth19_mod1    18 30351 30467 -15158    30315 129.89  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the 
# effect of GxE

# GxE is of higher magnitude effect when not accounting for family
anova(Growth19_mod1, Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## Growth19_mod1b   12 30469 30547 -15223    30445                         
## Growth19_mod1    18 30351 30467 -15158    30315 129.89  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT 
anova(Growth19_mod1,Growth19_mod1a,Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## Growth19_mod1a:     Family)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth19_mod1:     Family)
##                npar   AIC   BIC logLik deviance   Chisq Df Pr(>Chisq)    
## Growth19_mod1b   12 30469 30547 -15223    30445                          
## Growth19_mod1a   13 30358 30442 -15166    30332 113.392  1  < 2.2e-16 ***
## Growth19_mod1    18 30351 30467 -15158    30315  16.499  5   0.005555 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(Growth19_mod1,Growth19_mod1a,Growth19_mod1b)
##                df      AIC
## Growth19_mod1  18 30344.87
## Growth19_mod1a 13 30351.38
## Growth19_mod1b 12 30462.51
# effect size estimates for growth trait
plot_model(Growth19_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) + 
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

# diagnostic tests for the model
# sjPlot::plot_model(Growth19_mod1, type = "diag")

tab_model(Growth19_mod1)
  Growth
Predictors Estimates CI p
(Intercept) 9.53 7.86 – 11.20 <0.001
Garden [Maryland] 0.03 -2.24 – 2.31 0.976
Garden [North_Carolina] -0.41 -2.69 – 1.86 0.721
Region [Margin] 4.74 3.49 – 6.00 <0.001
Region [Edge] 3.02 2.05 – 3.98 <0.001
Garden [Maryland] *
Region [Margin] 		0.41 -0.79 – 1.61 0.506
Garden [North_Carolina] *
Region [Margin] 		-1.73 -2.95 – -0.51 0.005
Garden [Maryland] *
Region [Edge] 		0.29 -0.65 – 1.23 0.548
Garden [North_Carolina] *
Region [Edge] 		-0.96 -1.92 – -0.00 0.050
Random Effects
σ2 34.98
τ00 Family 3.55
τ00 Population 1.03
τ00 mBed 3.13
τ11 Family.GardenMaryland 0.89
τ11 Family.GardenNorth_Carolina 0.47
ρ01 Family.GardenMaryland 0.35
ρ01 Family.GardenNorth_Carolina -0.95
N Population 65
N mBed 15
N Family 334
Observations 4681
Marginal R2 / Conditional R2 0.089 / NA

Growth 2020

## working model for growth - ****
Growth_mod1 <- lmer(data = growth_2020, na.action = na.omit,
               Growth ~ Garden * Region +
               (1|Population) + (1|mBed) +
               (Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -2.3e+01
Growth_mod1a <- lmer(data = growth_2020, na.action = na.omit,
               Growth ~ Garden * Region +
               (1|Population) + (1|mBed) + (1|Family))

Growth_mod1b <- lmer(data = growth_2020, na.action = na.omit,
               Growth ~ Garden * Region +
               (1|Population) + (1|mBed))

# model 1 anova 
anova(Growth_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Garden        1602.55  801.27     2  13.77 17.6768 0.0001573 ***
## Region        2830.07 1415.03     2  67.79 31.2169 2.454e-10 ***
## Garden:Region  815.96  203.99     4 368.64  4.5002 0.0014624 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(Growth_mod1, Growth_mod1a) 
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## Growth_mod1a:     Family)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth_mod1:     Family)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)
## Growth_mod1a   13 30161 30244 -15067    30135                     
## Growth_mod1    18 30164 30280 -15064    30128 6.1533  5     0.2916
# effect of genotype
anova(Growth_mod1a, Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## Growth_mod1a:     Family)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## Growth_mod1b   12 30216 30293 -15096    30192                         
## Growth_mod1a   13 30161 30244 -15067    30135 57.164  1  4.009e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(Growth_mod1, Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth_mod1:     Family)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## Growth_mod1b   12 30216 30293 -15096    30192                         
## Growth_mod1    18 30164 30280 -15064    30128 63.317  6   9.51e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the 
# effect of GxE

# GxE is of higher magnitude effect when not accounting for family
anova(Growth_mod1, Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth_mod1:     Family)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## Growth_mod1b   12 30216 30293 -15096    30192                         
## Growth_mod1    18 30164 30280 -15064    30128 63.317  6   9.51e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT 
anova(Growth_mod1,Growth_mod1a,Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## Growth_mod1a:     Family)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## Growth_mod1:     Family)
##              npar   AIC   BIC logLik deviance   Chisq Df Pr(>Chisq)    
## Growth_mod1b   12 30216 30293 -15096    30192                          
## Growth_mod1a   13 30161 30244 -15067    30135 57.1641  1  4.009e-14 ***
## Growth_mod1    18 30164 30280 -15064    30128  6.1533  5     0.2916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(Growth_mod1,Growth_mod1a,Growth_mod1b)
##              df      AIC
## Growth_mod1  18 30157.95
## Growth_mod1a 13 30154.50
## Growth_mod1b 12 30209.53
# effect size estimates for growth trait
plot_model(Growth_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) + 
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

# diagnostic tests for the model
# sjPlot::plot_model(Growth19_mod1, type = "diag")

tab_model(Growth_mod1)
  Growth
Predictors Estimates CI p
(Intercept) 16.04 14.91 – 17.17 <0.001
Garden [Maryland] 0.87 -0.50 – 2.25 0.212
Garden [North_Carolina] -2.94 -4.33 – -1.56 <0.001
Region [Margin] 3.76 2.16 – 5.37 <0.001
Region [Edge] -0.16 -1.39 – 1.07 0.796
Garden [Maryland] *
Region [Margin] 		3.08 1.56 – 4.59 <0.001
Garden [North_Carolina] *
Region [Margin] 		2.69 1.12 – 4.26 0.001
Garden [Maryland] *
Region [Edge] 		0.49 -0.71 – 1.69 0.420
Garden [North_Carolina] *
Region [Edge] 		0.67 -0.57 – 1.92 0.289
Random Effects
σ2 45.33
τ00 Family 0.00
τ00 Population 3.10
τ00 mBed 0.86
τ11 Family.GardenMaryland 4.07
τ11 Family.GardenNorth_Carolina 6.30
ρ01 Family.GardenMaryland  
ρ01 Family.GardenNorth_Carolina  
N Population 65
N mBed 15
N Family 334
Observations 4475
Marginal R2 / Conditional R2 0.136 / NA

Bud Set 2019

budset_2019 <- budset_2019%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)

## working model for growth - ****
budset19_mod1 <- lmer(data = budset_2019, na.action = na.omit,
               BudSet ~ Garden * Region +
               (1|Population) + (1|mBed) +
               (Garden|Family)) # G x E
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Warning: Model failed to converge with 1 negative eigenvalue: -3.2e+01
budset19_mod1a <- lmer(data = budset_2019, na.action = na.omit,
               BudSet ~ Garden * Region +
               (1|Population) + (1|mBed) + (1|Family))

budset19_mod1b <- lmer(data = budset_2019, na.action = na.omit,
               BudSet ~ Garden * Region +
               (1|Population) + (1|mBed))

# model 1 anova 
anova(budset19_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Garden         39470 19735.0     2  13.60 31.1200 8.408e-06 ***
## Region         39318 19658.9     2  69.07 30.9999 2.465e-10 ***
## Garden:Region  12902  3225.6     4 604.80  5.0865 0.0004895 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budset19_mod1, budset19_mod1a) 
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budset19_mod1a:     Family)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset19_mod1a   13 44060 44144 -22017    44034                         
## budset19_mod1    18 44019 44135 -21991    43983 50.972  5  8.765e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budset19_mod1a, budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budset19_mod1a:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset19_mod1b   12 44135 44212 -22055    44111                         
## budset19_mod1a   13 44060 44144 -22017    44034 77.107  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budset19_mod1, budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset19_mod1b   12 44135 44212 -22055    44111                         
## budset19_mod1    18 44019 44135 -21991    43983 128.08  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the 
# effect of GxE

# GxE is of higher magnitude effect when not accounting for family
anova(budset19_mod1, budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset19_mod1b   12 44135 44212 -22055    44111                         
## budset19_mod1    18 44019 44135 -21991    43983 128.08  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT 
anova(budset19_mod1,budset19_mod1a,budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budset19_mod1a:     Family)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset19_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset19_mod1b   12 44135 44212 -22055    44111                         
## budset19_mod1a   13 44060 44144 -22017    44034 77.107  1  < 2.2e-16 ***
## budset19_mod1    18 44019 44135 -21991    43983 50.972  5  8.765e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budset19_mod1,budset19_mod1a,budset19_mod1b)
##                df      AIC
## budset19_mod1  18 43989.61
## budset19_mod1a 13 44029.72
## budset19_mod1b 12 44104.62
# effect size estimates for budset trait
plot_model(budset19_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) + 
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

tab_model(budset19_mod1)
  BudSet
Predictors Estimates CI p
(Intercept) 260.61 255.82 – 265.40 <0.001
Garden [Maryland] -8.61 -14.69 – -2.52 0.006
Garden [North_Carolina] -20.62 -26.83 – -14.40 <0.001
Region [Margin] 8.89 3.12 – 14.66 0.003
Region [Edge] 16.94 12.49 – 21.39 <0.001
Garden [Maryland] *
Region [Margin] 		-4.30 -9.64 – 1.05 0.115
Garden [North_Carolina] *
Region [Margin] 		-8.88 -15.02 – -2.74 0.005
Garden [Maryland] *
Region [Edge] 		-7.34 -11.55 – -3.13 0.001
Garden [North_Carolina] *
Region [Edge] 		-1.85 -6.66 – 2.96 0.450
Random Effects
σ2 634.16
τ00 Family 103.15
τ00 Population 20.11
τ00 mBed 19.75
τ11 Family.GardenMaryland 42.42
τ11 Family.GardenNorth_Carolina 100.27
ρ01 Family.GardenMaryland -1.00
ρ01 Family.GardenNorth_Carolina -0.39
ICC 0.16
N Population 65
N mBed 15
N Family 334
Observations 4685
Marginal R2 / Conditional R2 0.143 / 0.277

Bud Set 2020

budset_2020 <- budset_2020%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)


## working model for growth - ****
budset20_mod1 <- lmer(data = budset_2020, na.action = na.omit,
               BudSet ~ Garden * Region +
               (1|Population) + (1|mBed) +
               (Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
budset20_mod1a <- lmer(data = budset_2020, na.action = na.omit,
               BudSet ~ Garden * Region +
               (1|Population) + (1|mBed) + (1|Family))

budset20_mod1b <- lmer(data = budset_2020, na.action = na.omit,
               BudSet ~ Garden * Region +
               (1|Population) + (1|mBed))

# model 1 anova 
anova(budset20_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Garden         66924   33462     2  13.29 107.4274 6.252e-09 ***
## Region          6295    3148     2  62.25  10.1050 0.0001583 ***
## Garden:Region   3817     954     4 518.61   3.0635 0.0163764 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budset20_mod1, budset20_mod1a) 
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budset20_mod1a:     Family)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset20_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset20_mod1a   13 37141 37224 -18558    37115                         
## budset20_mod1    18 37116 37231 -18540    37080 35.178  5  1.386e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budset20_mod1a, budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budset20_mod1a:     Family)
##                npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)    
## budset20_mod1b   12 37185 37261 -18580    37161                        
## budset20_mod1a   13 37141 37224 -18558    37115 45.36  1   1.64e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budset20_mod1, budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset20_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset20_mod1b   12 37185 37261 -18580    37161                         
## budset20_mod1    18 37116 37231 -18540    37080 80.538  6  2.766e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the 
# effect of GxE

# GxE is of higher magnitude effect when not accounting for family
anova(budset20_mod1, budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset20_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset20_mod1b   12 37185 37261 -18580    37161                         
## budset20_mod1    18 37116 37231 -18540    37080 80.538  6  2.766e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT 
anova(budset20_mod1,budset20_mod1a,budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budset20_mod1a:     Family)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budset20_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budset20_mod1b   12 37185 37261 -18580    37161                         
## budset20_mod1a   13 37141 37224 -18558    37115 45.360  1  1.640e-11 ***
## budset20_mod1    18 37116 37231 -18540    37080 35.178  5  1.386e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budset20_mod1,budset20_mod1a,budset20_mod1b)
##                df      AIC
## budset20_mod1  18 37092.72
## budset20_mod1a 13 37118.40
## budset20_mod1b 12 37161.59
# effect size estimates for budset trait
plot_model(budset20_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) + 
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

tab_model(budset20_mod1)
  BudSet
Predictors Estimates CI p
(Intercept) 235.37 232.04 – 238.70 <0.001
Garden [Maryland] -25.29 -29.69 – -20.89 <0.001
Garden [North_Carolina] -28.81 -33.25 – -24.37 <0.001
Region [Margin] 6.71 2.89 – 10.53 0.001
Region [Edge] 2.02 -0.94 – 4.97 0.181
Garden [Maryland] *
Region [Margin] 		0.02 -4.04 – 4.08 0.991
Garden [North_Carolina] *
Region [Margin] 		-2.36 -6.50 – 1.79 0.266
Garden [Maryland] *
Region [Edge] 		-2.15 -5.34 – 1.04 0.187
Garden [North_Carolina] *
Region [Edge] 		2.00 -1.33 – 5.33 0.240
Random Effects
σ2 311.48
τ00 Family 57.16
τ00 Population 4.93
τ00 mBed 10.02
τ11 Family.GardenMaryland 37.27
τ11 Family.GardenNorth_Carolina 38.27
ρ01 Family.GardenMaryland -0.97
ρ01 Family.GardenNorth_Carolina -0.78
N Population 65
N mBed 15
N Family 334
Observations 4282
Marginal R2 / Conditional R2 0.361 / NA

Bud Break 2020

  • cGDD
budbreak_2020 <- budbreak_2020%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)

## working model for growth - ****
budbreak_mod1 <- lmer(data = budbreak_2020, na.action = na.omit,
               cGDD ~ Garden * Region +
               (1|Population) + (1|mBed) +
               (Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -6.2e-02 -6.7e+01
budbreak_mod1a <- lmer(data = budbreak_2020, na.action = na.omit,
               cGDD ~ Garden * Region +
               (1|Population) + (1|mBed) + (1|Family))

budbreak_mod1b <- lmer(data = budbreak_2020, na.action = na.omit,
               cGDD ~ Garden * Region +
               (1|Population) + (1|mBed))

# model 1 anova 
anova(budbreak_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Garden        2521553 1260777     2  17.95 1104.830 < 2.2e-16 ***
## Region         122334   61167     2  90.15   53.601 4.592e-16 ***
## Garden:Region   45697   11424     4 370.29   10.011 1.056e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budbreak_mod1, budbreak_mod1a) 
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budbreak_mod1a:     Family)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_mod1:     Family)
##                npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)    
## budbreak_mod1a   13 47668 47752 -23821    47642                        
## budbreak_mod1    18 47524 47641 -23744    47488 153.5  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budbreak_mod1a, budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budbreak_mod1a:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_mod1b   12 47825 47903 -23901    47801                         
## budbreak_mod1a   13 47668 47752 -23821    47642 159.29  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budbreak_mod1, budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_mod1b   12 47825 47903 -23901    47801                         
## budbreak_mod1    18 47524 47641 -23744    47488 312.79  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the 
# effect of GxE

# GxE is of higher magnitude effect when not accounting for family
anova(budbreak_mod1, budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_mod1b   12 47825 47903 -23901    47801                         
## budbreak_mod1    18 47524 47641 -23744    47488 312.79  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT 
anova(budbreak_mod1,budbreak_mod1a,budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budbreak_mod1a:     Family)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_mod1:     Family)
##                npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_mod1b   12 47825 47903 -23901    47801                         
## budbreak_mod1a   13 47668 47752 -23821    47642 159.29  1  < 2.2e-16 ***
## budbreak_mod1    18 47524 47641 -23744    47488 153.50  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budbreak_mod1,budbreak_mod1a,budbreak_mod1b)
##                df      AIC
## budbreak_mod1  18 47491.26
## budbreak_mod1a 13 47635.36
## budbreak_mod1b 12 47792.33
# effect size estimates for budbreak trait
plot_model(budbreak_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) + 
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

tab_model(budbreak_mod1)
  cGDD
Predictors Estimates CI p
(Intercept) 562.01 557.99 – 566.04 <0.001
Garden [Maryland] 90.76 85.39 – 96.13 <0.001
Garden [North_Carolina] 120.18 113.90 – 126.45 <0.001
Region [Margin] 7.47 1.55 – 13.39 0.013
Region [Edge] 12.03 7.48 – 16.59 <0.001
Garden [Maryland] *
Region [Margin] 		-1.43 -8.77 – 5.90 0.702
Garden [North_Carolina] *
Region [Margin] 		-6.34 -16.28 – 3.60 0.211
Garden [Maryland] *
Region [Edge] 		10.39 4.63 – 16.15 <0.001
Garden [North_Carolina] *
Region [Edge] 		21.99 14.18 – 29.80 <0.001
Random Effects
σ2 1141.15
τ00 Family 0.00
τ00 Population 24.59
τ00 mBed 10.48
τ11 Family.GardenMaryland 102.53
τ11 Family.GardenNorth_Carolina 560.79
ρ01 Family.GardenMaryland  
ρ01 Family.GardenNorth_Carolina  
N Population 65
N mBed 15
N Family 334
Observations 4757
Marginal R2 / Conditional R2 0.725 / NA

Visaulization

Growth 2019 rxn

  • paper**
    Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
Growth19_mod_nMb <- lmer(data = growth_2019,
               Growth ~ Garden * Region +
               (1|Population) + 
               (Garden|Family))
## boundary (singular) fit: see ?isSingular
heightData19 <- growth_2019
heightData19$pred <- predict(Growth19_mod_nMb)

# merge climateNA data
heightData19 <- merge(x=heightData19,
                     y=Population_climateNA[,c("Population","PC1","MAT")],
                     by="Population")


# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G19 <- ggplot(heightData19,aes(Garden,Growth))+
    facet_grid(.~Region)
    # zero_margin

# Rxn - family
G19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
                  xlab("Garden site") + ylab("Growth in 2019 (in cm)") +
  scale_color_viridis(option = "A", direction = -1) +
  theme_bw() + theme(axis.text=element_text(size=14), 
                     axis.title=element_text(size=14,face="bold"),
                     legend.title = element_text(color = "black", size = 14),
                     legend.text = element_text(color = "black", size = 13))

Region plot

# Rxn - population
growth19_pop_rxn <- G19 + geom_line(aes(y=pred,group=Family,colour=Population)) +
                  xlab("Garden site") + ylab("Vertical height growth (in cm)") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))




## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(heightData19,
       aes(x=Garden, y=Growth, color=Region)) +
  stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) + 
  # Lines by species using grouping
  stat_summary(aes(group = Region), geom = "line", fun = mean) +
  ylab("Growth in 2019 (in cm)") + theme(legend.position="right") +
  theme_bw() + theme(axis.text=element_text(size=14), 
                       axis.title=element_text(size=14,face="bold"),
                       legend.title = element_text(color = "black", size = 14),
                       legend.text = element_text(color = "black", size = 13)) +
  scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

Growth 2020 rxn

  • paper**
    Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
Growth20_mod_nMb <- lmer(data = growth_2020,
               Growth ~ Garden * Region +
               (1|Population) + 
               (Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.2e+02
heightData20 <- growth_2020
heightData20$pred <- predict(Growth20_mod_nMb)

# merge climateNA data
heightData20 <- merge(x=heightData20,
                     y=Population_climateNA[,c("Population","PC1","MAT")],
                     by="Population")


# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G20 <- ggplot(heightData20,aes(Garden,Growth))+
    facet_grid(.~Region)
    # zero_margin

# Rxn - family
G20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
                  xlab("Garden site") + ylab("Growth in 2020 (in cm)") +
  scale_color_viridis(option = "A", direction = -1) +
  theme_bw() + theme(axis.text=element_text(size=14), 
                     axis.title=element_text(size=14,face="bold"),
                     legend.title = element_text(color = "black", size = 14),
                     legend.text = element_text(color = "black", size = 13))

Region plot

# Rxn - population
growth20_pop_rxn <- G20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
                  xlab("Garden site") + ylab("Vertical height growth (in cm)") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))




## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(heightData20,
       aes(x=Garden, y=Growth, color=Region)) +
  stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) + 
  # Lines by species using grouping
  stat_summary(aes(group = Region), geom = "line", fun = mean) +
  ylab("Growth in 2020 (in cm)") + theme(legend.position="right") +
  theme_bw() + theme(axis.text=element_text(size=14), 
                       axis.title=element_text(size=14,face="bold"),
                       legend.title = element_text(color = "black", size = 14),
                       legend.text = element_text(color = "black", size = 13)) +
  scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

Budset 2019

Family plot

### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budset19_mod_nMb <- lmer(data = budset_2019,
               BudSet ~ Garden * Region +
               (1|Population) + 
               (Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -8.5e+00
budset_2019_2 <- budset_2019
budset_2019_2$pred <- predict(budset19_mod_nMb)

# merge climateNA data
budset_2019_2 <- merge(x=budset_2019_2,
                     y=Population_climateNA[,c("Population","PC1","MAT")],
                     by="Population")


# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BS19 <- ggplot(budset_2019_2,aes(Garden,BudSet))+
    facet_grid(.~Region)
    # zero_margin

# Rxn - family
BS19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
                  xlab("Garden site") + ylab("Bud set in 2019 (DOY)") +
  scale_color_viridis(option = "A", direction = -1) +
  theme_bw() + theme(axis.text=element_text(size=14), 
                     axis.title=element_text(size=14,face="bold"),
                     legend.title = element_text(color = "black", size = 14),
                     legend.text = element_text(color = "black", size = 13))

Region plot

# Rxn - population
budset19_pop_rxn <- BS19 + geom_line(aes(y=pred,group=Family,colour=Population)) +
                  xlab("Garden site") + ylab("Bud set in 2019 (DOY)") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))


## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(budset_2019_2,
       aes(x=Garden, y=BudSet, color=Region)) +
  stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) + 
  # Lines by species using grouping
  stat_summary(aes(group = Region), geom = "line", fun = mean) +
  ylab("Bud set in 2019 (DOY)") + theme(legend.position="right") +
  theme_bw() + theme(axis.text=element_text(size=14), 
                       axis.title=element_text(size=14,face="bold"),
                       legend.title = element_text(color = "black", size = 14),
                       legend.text = element_text(color = "black", size = 13)) +
  scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

Budset 2020

Family plot

### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budset20_mod_nMb <- lmer(data = budset_2020,
               BudSet ~ Garden * Region +
               (1|Population) + 
               (Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -6.8e+01
budset_2020_2 <- budset_2020
budset_2020_2$pred <- predict(budset20_mod_nMb)

# merge climateNA data
budset_2020_2 <- merge(x=budset_2020_2,
                     y=Population_climateNA[,c("Population","PC1","MAT")],
                     by="Population")


# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BS20 <- ggplot(budset_2020_2,aes(Garden,BudSet))+
    facet_grid(.~Region)
    # zero_margin

# Rxn - family
BS20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
                  xlab("Garden site") + ylab("Bud set in 2020 (DOY)") +
  scale_color_viridis(option = "A", direction = -1) +
  theme_bw() + theme(axis.text=element_text(size=14), 
                     axis.title=element_text(size=14,face="bold"),
                     legend.title = element_text(color = "black", size = 14),
                     legend.text = element_text(color = "black", size = 13))

Region plot

# Rxn - population
budset20_pop_rxn <- BS20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
                  xlab("Garden site") + ylab("Bud set in 2020 (DOY)") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))


## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(budset_2020_2,
       aes(x=Garden, y=BudSet, color=Region)) +
  stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) + 
  # Lines by species using grouping
  stat_summary(aes(group = Region), geom = "line", fun = mean) +
  ylab("Bud set in 2020 (DOY)") + theme(legend.position="right") +
  theme_bw() + theme(axis.text=element_text(size=14), 
                       axis.title=element_text(size=14,face="bold"),
                       legend.title = element_text(color = "black", size = 14),
                       legend.text = element_text(color = "black", size = 13)) +
  scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

Bud Break 2020

Family plot

### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budbreak_mod_nMb <- lmer(data = budbreak_2020,
               cGDD ~ Garden * Region +
               (1|Population) + 
               (Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.3e+02
budbreak <- budbreak_2020[!is.na(budbreak_2020$cGDD),]
budbreak$pred <- predict(budbreak_mod_nMb)

# merge climateNA data
budbreak <- merge(x=budbreak,
                   y=Population_climateNA[,c("Population","PC1","MAT")],
                   by="Population")

# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BB20 <- ggplot(budbreak,aes(Garden,cGDD))+
    facet_grid(.~Region)
    # zero_margin

# Rxn - family
BB20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
                  xlab("Garden site") + ylab("Bud break in 2020 (cGDD)") +
  scale_color_viridis(option = "A", direction = -1) +
  theme_bw() + theme(axis.text=element_text(size=14), 
                     axis.title=element_text(size=14,face="bold"),
                     legend.title = element_text(color = "black", size = 14),
                     legend.text = element_text(color = "black", size = 13))

Region plot

# Rxn - population
budbreak_pop_rxn <- BB20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
                  xlab("Garden site") + ylab("cGDD") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))


## Reaction norm for geo genetic groups/regions
ggplot(budbreak,
       aes(x=Garden, y=cGDD, color=Region)) +
  stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) + 
  # Lines by species using grouping
  stat_summary(aes(group = Region), geom = "line", fun = mean) +
  ylab("Bud break in 2020 (cGDD)") + theme(legend.position="right") +
  theme_bw() + theme(axis.text=element_text(size=14), 
                       axis.title=element_text(size=14,face="bold"),
                       legend.title = element_text(color = "black", size = 14),
                       legend.text = element_text(color = "black", size = 13)) +
  scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

Bud Break 2020 - cCDD

Models

budbreak_cCDD <- read.csv("C:\\Dropbox\\UVM\\Research\\Code/CommonGardenTraits/2020/BudBreak_2020_cCDD.csv", header=T,
                        stringsAsFactors = T)
# reordering factors
budbreak_cCDD$Region <- factor(budbreak_cCDD$Region,
                       levels = c("Core", "Margin", "Edge"))
budbreak_cCDD$Garden <- factor(budbreak_cCDD$Garden,
                       levels = c("Vermont", "Maryland", "North_Carolina"))
budbreak_cCDD <- budbreak_cCDD%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)





## working model for cCDD - ****
budbreak_cCDD_mod1 <- lmer(data = budbreak_cCDD,
               cCDD ~ Garden * Region +
               (1|Population) + (1|mBed) +
               (Garden|Family), na.action = na.omit) # G x E
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -8.8e-02 -9.2e-01
budbreak_cCDD_mod1a <- lmer(data = budbreak_cCDD,
               cCDD ~ Garden * Region +
               (1|Population) + (1|mBed) + (1|Family), na.action = na.omit)

budbreak_cCDD_mod1b <- lmer(data = budbreak_cCDD,
               cCDD ~ Garden * Region +
               (1|Population) + (1|mBed), na.action = na.omit)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.09547 (tol = 0.002, component 1)
# model 1 anova 
anova(budbreak_cCDD_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF    F value    Pr(>F)    
## Garden        125141   62570     2  21.26 83756.3158 < 2.2e-16 ***
## Region            34      17     2 400.59    22.5735 5.136e-10 ***
## Garden:Region     29       7     4 485.70     9.5912 1.801e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budbreak_cCDD_mod1, budbreak_cCDD_mod1a) 
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1a: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budbreak_cCDD_mod1a:     Family)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_cCDD_mod1:     Family)
##                     npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_cCDD_mod1a   13 13162 13246 -6568.0    13136                         
## budbreak_cCDD_mod1    18 12728 12845 -6346.3    12692 443.39  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budbreak_cCDD_mod1a, budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1a: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budbreak_cCDD_mod1a:     Family)
##                     npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_cCDD_mod1b   12 13230 13308 -6603.1    13206                         
## budbreak_cCDD_mod1a   13 13162 13246 -6568.0    13136 70.354  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budbreak_cCDD_mod1, budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_cCDD_mod1:     Family)
##                     npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_cCDD_mod1b   12 13230 13308 -6603.1    13206                         
## budbreak_cCDD_mod1    18 12728 12845 -6346.3    12692 513.74  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the 
# effect of GxE

# GxE is of higher magnitude effect when not accounting for family
anova(budbreak_cCDD_mod1, budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_cCDD_mod1:     Family)
##                     npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## budbreak_cCDD_mod1b   12 13230 13308 -6603.1    13206                         
## budbreak_cCDD_mod1    18 12728 12845 -6346.3    12692 513.74  6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT 
anova(budbreak_cCDD_mod1,budbreak_cCDD_mod1a,budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1a: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | 
## budbreak_cCDD_mod1a:     Family)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | 
## budbreak_cCDD_mod1:     Family)
##                     npar   AIC   BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## budbreak_cCDD_mod1b   12 13230 13308 -6603.1    13206                          
## budbreak_cCDD_mod1a   13 13162 13246 -6568.0    13136  70.354  1  < 2.2e-16 ***
## budbreak_cCDD_mod1    18 12728 12845 -6346.3    12692 443.389  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budbreak_cCDD_mod1,budbreak_cCDD_mod1a,budbreak_cCDD_mod1b)
##                     df      AIC
## budbreak_cCDD_mod1  18 12762.19
## budbreak_cCDD_mod1a 13 13196.03
## budbreak_cCDD_mod1b 12 13264.19
# effect size estimates for budbreak_cCDD trait
plot_model(budbreak_cCDD_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) + 
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

Reaction norm - part 1

# reaction norms for interactions
sjPlot::plot_model(budbreak_cCDD_mod1, type = "int")

sjPlot::plot_model(budbreak_cCDD_mod1, type = "re")
## [[1]]

## 
## [[2]]

## 
## [[3]]

##### reaction norms based on the model - Family
budbreak_cCDD3 <- budbreak_cCDD[!is.na(budbreak_cCDD$cCDD),]

# budbreak_cCDD2$pred <- predict(budbreak_cCDD_mod1, re.form=NULL) # all random effects are included

budbreak_cCDD3$pred <- predict(budbreak_cCDD_mod1,re.form=NA)  ## population level
budbreak_cCDD3$pred1 <- predict(budbreak_cCDD_mod1) ## individual level

g0 <- ggplot(budbreak_cCDD3,aes(Garden,cCDD, label=Family))+
    geom_point()+
    geom_line(aes(group=Family)) +
    facet_grid(.~Region)
    # zero_margin
g0 +   geom_line(colour="blue",aes(y=pred1,group=Family)) +
    geom_line(colour="red",aes(y=pred,group=Family)) 

##### reaction norms based on the model - Population

g0 <- ggplot(budbreak_cCDD3,aes(Garden,cCDD))+
    geom_point()+
    geom_line(aes(group=Population)) +
    facet_grid(.~Region)
    # zero_margin
g0 +   geom_line(colour="blue",aes(y=pred1,group=Population)) +
    geom_line(colour="red",aes(y=pred,group=Population))

Reaction norm - Part 2

# reaction norms based on the model - Family
budbreak_cCDD3 <- budbreak_cCDD[!is.na(budbreak_cCDD$cCDD),]

# budbreak_cCDD2$pred <- predict(budbreak_cCDD_mod1, re.form=NULL) # all random effects are included

budbreak_cCDD3$pred <- predict(budbreak_cCDD_mod1,re.form=NA)  ## population level/ no random effects included
budbreak_cCDD3$pred1 <- predict(budbreak_cCDD_mod1) ## individual level

G0 <- ggplot(budbreak_cCDD3,aes(Garden,cCDD))+
    # geom_point()+
    # geom_line(aes(group=Family)) +
    facet_grid(.~Region) #+ zero_margin


G0 + geom_line(colour="red",aes(y=pred,group=Family))  +
                  xlab("Garden site") + ylab("Vertical height growth (in cm)") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

budbreak_cCDD 2020 rxn - paper

### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budbreak_cCDD_mod_nMb <- lmer(data = budbreak_cCDD,
               cCDD ~ Garden * Region +
               (1|Population) + 
               (Garden|Family), na.action = na.omit)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -1.0e+00 -1.1e+00
budbreak_cCDD6 <- budbreak_cCDD3
budbreak_cCDD6$pred <- predict(budbreak_cCDD_mod_nMb)

# merge climateNA data
budbreak_cCDD6 <- merge(x=budbreak_cCDD6,
                   y=Population_climateNA[,c("Population","PC1","MAT")],
                   by="Population")

# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G0 <- ggplot(budbreak_cCDD6,aes(Garden,cCDD))+
    facet_grid(.~Region)
    # zero_margin

# Rxn - family
G0 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
                  xlab("Garden site") + ylab("Bud break in 2020 (cCDD)") +
  scale_color_viridis(option = "A", direction = -1) +
  theme_bw() + theme(axis.text=element_text(size=14), 
                     axis.title=element_text(size=14,face="bold"),
                     legend.title = element_text(color = "black", size = 14),
                     legend.text = element_text(color = "black", size = 13))

# Rxn - population
budbreak_cCDD_pop_rxn <- G0 + geom_line(aes(y=pred,group=Family,colour=Population)) +
                  xlab("Garden site") + ylab("cCDD") +
  theme_minimal() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"),
        strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))

plotly::ggplotly(budbreak_cCDD_pop_rxn)
budbreak_cCDD_pop_rxn

## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(budbreak_cCDD3, aes(x = Garden, y = cCDD, color = Region)) +
  stat_summary(aes(group = Region), fun = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun = mean, geom = "point", size = 4) +
  # geom_point(aes(color = Region), size = 2) +
  theme_minimal() + ylab("cCDD") + 
  theme(text = element_text(size=12),
        axis.text.x = element_text(angle=0))

Reaction norm - Part 3

ggplot(budbreak_cCDD3,
       aes(x=Garden, y=cCDD, color=Region)) +
  stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
  stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
  stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) + 
  # Lines by species using grouping
  stat_summary(aes(group = Region), geom = "line", fun = mean) +
  ylab("Bud break in 2020 (cCDD)") + theme(legend.position="right") +
  theme_bw() + theme(axis.text=element_text(size=14), 
                       axis.title=element_text(size=14,face="bold"),
                       legend.title = element_text(color = "black", size = 14),
                       legend.text = element_text(color = "black", size = 13)) +
  scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

# diagnostic tests for the model
sjPlot::plot_model(budbreak_cCDD_mod1, type = "diag")
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$Family
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]$Population
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]$mBed
## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

tab_model(budbreak_cCDD_mod1)
  cCDD
Predictors Estimates CI p
(Intercept) 80.26 80.17 – 80.35 <0.001
Garden [Maryland] 13.54 13.37 – 13.71 <0.001
Garden [North_Carolina] -14.78 -14.91 – -14.64 <0.001
Region [Margin] -0.01 -0.13 – 0.11 0.875
Region [Edge] -0.00 -0.10 – 0.09 0.951
Garden [Maryland] *
Region [Margin] 		-0.06 -0.35 – 0.23 0.698
Garden [North_Carolina] *
Region [Margin] 		-0.01 -0.20 – 0.17 0.888
Garden [Maryland] *
Region [Edge] 		0.49 0.26 – 0.71 <0.001
Garden [North_Carolina] *
Region [Edge] 		0.42 0.27 – 0.57 <0.001
Random Effects
σ2 0.75
τ00 Family 0.00
τ00 Population 0.00
τ00 mBed 0.01
τ11 Family.GardenMaryland 0.57
τ11 Family.GardenNorth_Carolina 0.06
ρ01 Family.GardenMaryland  
ρ01 Family.GardenNorth_Carolina  
N Population 65
N mBed 15
N Family 334
Observations 4757
Marginal R2 / Conditional R2 0.994 / NA